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Abstract 

This paper studies and analyzes a preconditioned Krylov solver for Helmholtz problems that are 
formulated with absorbing boundary layers based on complex coordinate stretching. The preconditioner 
problem is a Helmholtz problem where not only the coordinates in the absorbing layer have an imaginary 
part, but also the coordinates in the interior region. This results into a preconditioner problem that is 
invertible with a multigrid cycle. We give a numerical analysis based on the eigenvalues and evaluate the 
performance with several numerical experiments. The method is an alternative to the complex shifted 
Laplacian and it gives a comparable performance for the studied model problems. 

1 Introduction 

The Helmholtz equation is frequently used to model propagation of waves in applications such as acoustic, 
seismic and electromagnetic realistic systems. It is less known that the equation is also helpful to understand 
and predict the reaction rates of fundamental processes in few-body physics and chemistry that are important 
for many areas of technology. In gas discharge reactors that are used industrially for chemical processing 
of surfaces, for example, these reaction rates are an essential modeling input [1]. There is both a scientific 
interest and an industrial need for accurate prediction of reaction rates [2] . To predict accurately the reaction 
rates of these processes it is necessary to solve the multi-dimensional Schrodinger equation [3, 4] that is, for 
the energy regime of these processes, equivalent to a multi-dimensional Helmholtz equation with outgoing 
waves boundary conditions and a space dependent wave number. The reaction rate of a particular process 
is then found as a post-processing step where the fluxes of the outgoing waves, corresponding to a particular 
reaction, are extracted [5]. 

These applications, often in diverse fields, have little in common amongst them except the Helmholtz model 
used in the simulations. Naturally, the numerical solution of the indefinite Helmholtz equation forms an 
interesting field of research for a widespread scientific community. Two main numerical concerns in this 
context are the truncation of the infinite physical domain to a finite numerical one mapped on a grid, and 
the efficient iterative solution of the resulting indefinite discrete linear system. In this paper, we concentrate 
on the latter of these issues. 

A hard truncation of the physical domain (in the Dirichlet sense) on a numerically feasible finite boundary, 
results in waves reflecting back and propagating into the truncated domain. These artificial reflections 
have to be avoided for an acceptable numerical treatment of the Helmholtz equation. This consideration 
led to the commonly used boundary conditions by Engquist and Majda [6] and Bayliss and Turkel [7] 
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who approximated the Sommerfeld radiation condition for homogeneous media at the boundary. Berenger 
avoided reflecting waves by adding perfectly matched layers (PML) [8] to the truncated domain. On this 
absorbing extension the problem is reformulated in order to damp the solution exponentially. Both techniques 
have been generalized and fine-tuned by many authors, see e.g. [9] for an overview. Chew and Weedon 
[10] related the PML method to a complex coordinate stretching. They consider the absorbing layer as 
an analytical continuation of the space domain into the complex plane, where the original equation is 
preserved. Earlier, in the 70's, linear complex scaling of the space domain was already used as a method 
to compute atomic resonances in microscopic system described by the Schrodinger equation [11, 12]. There 
a coordinate transformation, r — >• rexp(z#) with a angle > 0, was applied on the full domain. In the 
same decade, Simon introduced exterior complex scaling (ECS) by only transforming the boundary region 
with the transformation r — ^ (Rq — r) exp(i#) + r, where Ro denotes the start of the boundary region. The 
purpose of this transformation was to introduce the correct boundary condition for resonant states while 
avoiding analytical continuation of non-analytical potentials in the Schrodinger equation. Later this ECS 
transformation was used to enforce outgoing wave boundary conditions to atomic break-up problems [3]. 
Note that the ECS transformation has a discontinuous first derivative in i?o? while in PML the complex 
stretching is usually introduced as a smooth coordinate transformation. 

The purpose of this paper is to find efficient iterative solvers for Helmholtz problems equipped with absorbing 
boundary layers based on complex coordinate stretching. The main iterative challenge in a discrete Helmholtz 
problem H^Uh = bh is the indefiniteness of the discretized operator H^. A powerful way to get around 
this issue is the use of preconditioned Krylov subspace methods. The original troublesome matrix is 
multiplied by the inverse of the preconditioning matrix M^, resulting in a new system M^Hh = M^ 1 bh 
for left preconditioning. The choice of the preconditioner is a trade off between a cheaply invertible matrix 
Mh and a definite preconditioned system M^H^ with nicely clustered eigenvalues. The former demand 
can be relaxed by allowing an inexact inversion of the preconditioner M^, e.g. a few sweeps of a multigrid 
method. A successful preconditioner in this setup of multigrid preconditioning (MGP) of Krylov methods is 
the complex shifted Laplacian (CSL) developed for Sommerfeld radiation conditions by Erlangga, Vuik and 
Oosterlee [13]. 

In this paper we study the effect of complex stretching a part of the domain on the eigenvalues of the 
discretized Helmholtz equation. We have chosen to study the simplest problem with complex coordinate 
stretching which is the ECS domain as introduced by Simon and still used for atomic break-up problems. 
We have analyzed the eigenvalues of the one-dimensional Laplacian discretized on an ECS domain. Then the 
achieved insights are used to apply the CSL preconditioning idea to two-dimensional Helmholtz problems 
with ECS boundaries. The theoretical analysis also leads to an alternative family of preconditioners based on 
different complex stretchings of the numerical grid (CSG). Although ECS is not the most accurate absorbing 
boundary condition, we believe the insights on the performance of the iterative solver are valid for problems 
where a smooth complex stretching transformation is introduced. 

The remainder of this paper is organized as follows. In Section 2, we describe the ECS absorbing boundary 
layer. We use a one-dimensional Laplace problem for theoretical considerations; this reference problem is also 
described in Section 2. The numerical analysis of the discrete one-dimensional reference problem is given next 
in Section 3. Three important lemmas are given here, the insights from which paved the way for the work in 
this paper. We use three model problems for experimentation. They possess particular properties, which we 
detail in Section 4. Section 5 follows, and deals with preconditioning ideas. Multigrid behavior as a solver 
as well as for approximate preconditioner inversion is discussed. We also calibrate multigrid performance 
with different components. Numerical validation of all the theoretical insights is given in Section 6 where 
the model problems are solved with multigrid preconditioned Bi-CGSTAB [14] and IDR(«s) [15] (with s = 4 
and 5 = 8), with the CSL and the CSG preconditioning operators, and accompanied by comparison tables. 
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2 Exterior complex stretched domains 



The Helmholtz equation in a homogeneous medium is 
Definition 2.1 (Helmholtz). 

Hu = - (A + k 2 ) u = x in Q K d , (2.1) 

with dimension d > 1. 

where x is a source term, k G R is called the wave number. The equation describes acoustic wave problems 
on an unbounded domain Qq. We want to solve the Helmholtz equation (2.1) numerically on a bounded part 
of the domain Q C ^o, with absorbing boundary layers. 

For the Helmholtz equation (2.1) restricted to the bounded domain £7, a popular boundary condition is the 
first order Sommerfeld boundary condition, given by, 

Ou 

—— — —iku on dQ, (2.2) 
on 

where dft represents the domain boundary and n, the outward normal. We write i for the complex identity. In 
multi-dimensional Helmholtz problems absorbing boundary layers are preferred over these classical first order 
Sommerfeld conditions because the latter requires an exact knowledge of the wave number at the boundary. 
More important, in higher dimensions condition (2.2) suffers from artificial reflections and a higher order 
version should be applied [6, 7]. Equation (2.2) is still useful for the analysis of iterative methods though. 
The physical interpretation of absorbing boundary layers is to extend the original domain Q with a layer T 
of an absorbing material. In Q the original equation is kept and in the layer T the equation is manipulated 
to enforce specific boundary conditions on the new boundaries ffT, through an adapted potential due to 
a change in the material. This was introduced as a perfectly matched layer (PML) by Berenger [8]. The 
original PML idea is mathematically equivalent to a particular complex coordinate stretching [10] in the 
boundary layers, where the original equation is used in a new coordinate system. In this complex stretching 
approach we define an analytic continuation on the layers by 

I x + ij(x), x G 1 , 

with / G C 2 the stretching function, increasing (e.g. linear, quadratic, . . . ) and lim f(x) = 0. We denote 
the image of the layer T z = z(T) and call it the complex contour. These robust boundary layers do not 
use the wave number explicitly as opposed to the Sommerfeld conditions (2.2) and they can easily be tuned 
in numerical experiments. The transformation (2.3) constructs the absorbing complex contour by adding 
a complex shift to the domain extension T, but it can also be done by a complex rotation of T. For a 
one-dimensional problem in Q = [xo,r] with linear complex stretching applied to the extension T — [r, R], 
with an angle this is typically 

z(x) = (x — r)e %e + r, x G T, 

and was introduced by Simon as exterior complex scaling [16]. On discrete level the mesh width on the 
contour T z becomes h 1 = he %d . We point out that although the above definition suffices well for practical 
purposes, as for the experiments presented in this paper, we will keep the analytic discussion general by 
using the expression (2.3). 

We will study the effect of a complex stretching transformation on a simple one-dimensional Laplace problem 

d 2 

- Lu(x) = --j^u(x) = x(x) in [0, 1] C R, (2.4) 

with u(0) = and an absorbing boundary condition in x = 1. The minus sign is introduced to make the 
Laplacian positive definite. The homogeneous Helmholtz problem only differs in a constant shift k 2 and 
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shares the same eigenvectors. An eigenvalue of the Helmholtz problem is Xl — k 2 where A^ is an eigenvalue 
of the negative Laplacian L. However, this minor difference can turn the definite problem into a indefinite 
one, with major consequences on the behavior of iterative methods such as Krylov methods and multigrid 
methods. The numerical results in Section 6 show that the results are easily extended to two dimensions, 
and that the same preconditioning ideas can be useful for different values of k and for non- homogeneous 
Helmholtz problems. 

We now implement the complex stretched boundary layer on the one-dimensional Laplace problem (2.4) by 
adding an extension T = [1,R] with 1 < R G M to construct the complex contour T z = z(T) C C. In this 
paper, we use a linear coordinate transformation on the layer so that Y z is the complex line connecting 1 
and z(R) = R z G C that we will denote as = [l,R z ). This transforms the problem to 

- Lu(z) = -^u(z) = x (z) in [0, 1] U [1, R z ) C C, (2.5) 

with homogeneous Dirichlet boundary conditions at z(0) = and z(R) = R z (see Figure 1). In the re- 
mainder of the paper will refer to this linear stretching (2.5) as the exterior complex scaled (or stretched) 
transformation or in short ECS. The two boundary points z = and z = R z determine the position of the 
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Figure 1: The ECS domain z(x). An ECS contour is added as an extension of the domain at x = 1. This 
shifts the domain and consequently the spectrum of the resulting operator into the complex plane. 
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eigenvalues independently of the path connecting the two points. This is readily obvious by inspecting the 
following equation, which gives the eigenvalues of the negative Laplacian on any one-dimensional curve in 
the complex plane connecting the points and R z : 

Al= (/?) with ^' GN o- 

The derivation is trivial and hence not shown here. The complex contour acts as an absorbing layer. Indeed, 
because of the exponential decay of the analytically continued solution one can enforce homogeneous Dirichlet 
boundary conditions at the end of the complex contour with a significantly smaller boundary error than on 
a real truncated domain [17], as is illustrated in Figure 2. The larger the imaginary part of the complex 
boundary, the stronger the suppression of the reflected waves. Although the shape of the contour, i.e. the 
stretching function f(x), does not explicitly influence the damping in the continuous case, the discretized 
problem is susceptible to different shapes, as we will see in the next section. 

Remark 1 (Accuracy of ECS layers). A detailed discussion on the accuracy of ECS layers does not lie 
within the scope of this paper; we merely use the simple linear ECS in (2.5) as a model to understand the 
iterative solution of the problems with more advanced complex stretched boundary layers. 



4 



real 



Figure 2: ECS on a one-dimensional domain. A complex contour is added to both boundaries, enforcing the 
wave to fulfill Dirichlet boundary conditions. e lKZ ^ = e m O+*/0)) = e iKx e -K,f(x) w ( co lor online) 



3 Numerical analysis of discretized operator 

In this section we discretize the Laplace problem (2.5) with finite differences using the Shortley-Weller 
formula for non-uniform vertex centered grids [18]. It enables discretization through the region of transition 
to complex mesh widths for the complex contour in the ECS domain. We present theoretical results for the 
eigenvalues of the discretization matrix. 

Consider the one-dimensional Laplace problem (2.5). We define a uniform grid 

(zj)o<j<n on [0, 1] 

with zo = and z n = 1 and mesh width h = 1/n G M, and a second uniform grid on the complex contour 

( z j)n<j<n+m On [l,i2^] 

with z n+m = R z and complex mesh width h 7 = (R z — l)/m. We will refer to the angle of h 7 in the complex 
plane as the ECS angle and denote it 9 1 . The union of these two grids is the ECS grid 

(%)o<j<n+m on [0, 1] U [1, R z ] (3.1) 

in the entire ECS domain. We will often use the fraction 7 = h 7 /h G C. To approximate the second 
derivative in (2.5) we choose the Shortley-Weller formula 

d 2 U/ s 2 ( \ 



for non-uniform grids in grid point j, where hj-i and hj are the left and right mesh widths respectively, 
and may belong either to the h category or to the h^ category. The formula is easily derived from Taylor 
series and reduces to regular second order central differences when hj-\ = hj, i.e., in the interior real region 
(0,1), and in the interior of the complex contour (1,R Z ) because the stretching function / is taken to be 
linear. The only exception is the point z n where at most we lose an order of accuracy, however with ample 
discretization steps, the overall accuracy is anticipated to match up to second order. The result is a linear 
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Figure 3: Discretized ECS domain Zj. The ECS domain is discretized with complex mesh widths on the 
complex contour. 



system of equations that we will represent by the matrix equation 



- L h u h = b h . 



(3.2) 



The right hand side bh contains contributions from the source function x- The spectrum of the discretization 
matrix Lh in (3.2) determines the convergence behavior of iterative methods such as Krylov subspace methods 
and multigrid schemes for solving the system. It is drastically different from the spectrum of the continuous 
operator L. We start with the construction of bounds on the field of values of — L^. 

The remainder of this section is focused on several lemmas that will help to understand the spectral properties 
of the discretized operator. First, in Lemma 3.1 the Gershgorin disks are used to produce bounds on the 
spectrum. Next, in Lemma 3.2 we find a condition for the eigenvalues of the discrete ECS Helmholtz operator 
for constant wave numbers. The solutions lie on a pitchfork-shaped figure. In Lemma 3.3 of this section, it 
is shown how an approximation can be found for the limiting points of this spectrum. 

Lemma 3.1. Consider the ECS grid (3.1) and the discretization matrix Lh in Equation (3.2). Define 
7 = and its complex conjugate 7. If A £ cr(-Lh), then 
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< h 2 $i(\) <max(4,3 



3 + 7 



1 + 7 



T 



< 



|7| 5 



l7| 2 +7 



where 5ft and 9 denote the real and imaginary part, respectively. 

Proof. Every eigenvalue of a matrix lies in the field of values of that matrix. For the field of values W(A^) 
of a matrix Am = {aij)i<i,j<N holds 

mm{^{W{A N ))} = min{fi€R: fi€a( AN ^ A * N )} 

max{St(W(A N ))} = max{/x g R : fi € <r( — ± A * N )} 

mm{3(W{A N ))} = min{fi eR:[i£ -ur( ~ N )} 

max{Q(W(A N ))} = maxfjit G R : fi G -kj{ An ~ An )}, 



with A* N the adjoint matrix of Am- We scale the model problem with ft, 2 and construct the Gershgorin disks 
of —ft 2 Lh l" Lfe and —ft 2 Lh Z Lh ■ The eigenvalues of a matrix lie on the union of its Gershgorin disks. They 
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are defined as 



The scaled operator is 



N 



Vi = {zeC: \z-CLjjl <J2\ a ij\}' 
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We will use the notation B(c,r) C C for the ball centered around c G C with radius r G 
-h 2 ^^ has seven distinct Gershgorin disks 



i = 1 : V 1 
M2<i<n-2:V 2 
i = n — 1 : V3 

i = n : V4 

j = m - 1 : £> 5 
V2<j<ra-2:£> 6 
j = 1 : D 7 



5(2,1) 
5(2,2) 



1+7 



B 25ft( T 
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So the matrix 



The minimum and maximum of (J 1<i<7 = U2<i<6^* determine a lower and upper bound for h 2 cr( 
For our ECS problems, ECS angles up to | (0 < 6 1 < |) the minimum is J?(^2 ) - 
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maximum is max (4, 3+1 
and does not exclude negative values 
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). Note that the lower bound for the real part of the eigenvalues is negative 



The matrix —h 
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has five distinct Gershgorin disks 
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The minimum and maximum of Ui<i<5 ^ = Ui<i<4 ^ determine a lower and upper bound for —ih 2 a(— Lh 2 Lh ). 
For our ECS problems, ECS angles up to | (0 < 6 1 < |) these extrema are 4^(^2 ) and | 
tively. □ 
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Figure 4: Spectrum of — h 2 Lh (•) and the Gershgorin disks of —h 2 /2(Lh + (red dotted circles) and 
—h 2 /2 (Lh — L* h ) (blue dotted circles). The spectrum is bounded by the field of values that lies inside the 
rectangle derived from the Gershgorin disks. The red circles lead to the left and right bound, the blue circles 
give the upper and lower bound. The negative left bound does not assure positive definiteness of the matrix, 
(color online) 



By considering the Gershgorin disks of only the discretized Laplacian these bounds can be further 
sharpened. Next is a discussion on the exact position of the eigenvalues. 

Lemma 3.2. Consider the ECS grid (3.1) and the discretization matrix in Equation (3.2). Define 
7=7^- Then the eigenvalues of —L^ are the solutions of 

tan(2np(A)) cosj^ = 
v ; tan(2rag(A)) cos(<?(A)) ' v ; 

with p(X) = \ arccos(l — ^h 2 ), q(X) — \ arccos(l — ^j 2 h 2 ). 

Proof Note that the eigenvalues of the Helmholtz problem fulfill the same condition (3.3), with \n h = 
\L h — k 2 . We prove the Laplace case. 

We write the grid points slightly different than before by numbering the grid points {xj)\< n in [0, 1] from 
left to right, and (yj)i<j<m in [1?^] form right to left, so that the turning point x n = y m = 1. In other 
words we consider the ECS grid as two joint grids (xj)i<j<n an d (%)i<j<m- Consequently the two Dirichlet 
boundaries are xo and yo. The Shortley-Weller finite differences formula reduces to regular second order 
central differences in every grid point, except for the turning point x n = y m = 1 with left grid distance h and 
right grid distance jh. We look for an eigenvector v on the real grid and an eigenvector w on the complex 



8 



contour with the same eigenvalue A and v n = w m in the turning point. We get the recurrence relations 

-^2 (vj-i - 2vj + v j+1 ) = Xvj, < j < n; 

(wj-i - 2^- + Wj+i) = A^-, < j < m. 

The first and the last are Chebyshev recurrence relations with general solutions 

j vj = ciVj(sx) + c 2 Tj{s\), < j < n; 
[ ^ = hVjitx) + d 2 Tj(t x ), < j < m. 

with 8a = 1 — f = 1 — f (7^) 2 an d V}(sa) = VI — s \Uj-i(s\). Tj and J7j are the j-th Chebyshev 

polynomials of the first and the second kind respectively. We apply the boundary conditions uo = and 
= and find 

f Vj = C!Vj(s\), < j < n; 
[ ^ = diV^A), < j < m. 

The matching condition returns 

Vn = 

V»(s A ) 



where we assumed v n 7^ so V n (s \) / / V m (t\). The remaining constant ci determines the norm of 
the eigenvector. So we can assume c\ = 1. Now the only free parameter left is the eigenvalue A that is 
determined by the second recurrence relation, the discretization scheme in the turning point. 

Vn-M - ( 1 + -) V n (s x ) + Ig^4^ m _i(t A ) > ) = AK( SA ) 



^ 2 (l + 7) V V 7 V 77 V AJ lV m (t x ) 

^(l + 7 ) V K( SA ) V 7/ 7 V m (t x ) J 

^ Vm-litx) , K-i(*a) 
Kn(^A) ^ n (s A ) 

sin((ra — 1) arccos(tA)) sin((n — 1) arccos(sA)) 

& —( ~u\\ h ^ ~t ? — v> — = 7 5 a+^a 

sm m arccos(tAj) sin n arccos(sAj) 



t\ - cot(2mq x )\ 1 - t\ + 7 [ s x - cot(2np x )\ 1 - s% = 7^ A + t A 



cot(2mqx)\/ 1 — t x — J cot(2npx) \ 1 — s 



tan(2npA) / 1 — 5 



tan(2rag A ) V 1 - £ 2 



A 



tan(2np A ) /1 + s A 



tan(2m<7A) V 1 + ^a 



We introduced the shorthands p x = arccos(t) and q x = arccos(s), substituted Vj(x) = sin(j arccos(x)) and 
excluded the trivial cases cot(2npA) = and t\ = l. □ 

We can now solve the eigenvalue problem numerically by applying e.g. Newton's method on the function 
(3.3). There are eigenvalues to be found along the complex line pe~ 2z6 ^ with p G R, and close to A/h 2 and 
4/ 7 2 /i 2 . 
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Lemma 3.3. Let —L^ be the negative discretized Laplacian in Equation (3.2) with eigenvalues A G o~(— Lh). 
Then three typical regions of the eigenvalues can be identified in the spectrum: 

4 7T 

For |A - — j I < 1 : A w 4n 2 sin(Z — ) tu^ft Z < n 
Zr 2n 

4 7T 

i^or I A - — ^ — ^ I < 1 : A w 47 2 m 2 sin(Z ) wit/i Z < m 

For |A| « 1 : A « f //•//// / > 1 



These approximations are the largest eigenvalues of the discretized Laplacian restricted to the real domain 
[0,1], the complex contour [1,R Z ], and the smallest eigenvalues of the discretized Laplacian on the complex 
line [Q,R Z ], respectively. 

Proof We take the scaled operator —h 2 Lh- The eigenvalues are the roots of the function 

F(fi) = sin (2np(/i)) cos (2mg(/i)) cos (#(//)) (3.4) 
+ cos (2np(/i)) sin (2mg(/x)) cos (p(/J>)) 

with = | arccos(l — §), = \ arccos(l — ^7 2 ). Using Taylor series we get 

P(M) = |-^V / 4^+0(|4-M| 3/2 ) 
for ji « 4, so .F(/i) can be approximated by 

F(/i) « sin ^7r — ^4 — cos (2mq{n)) cos (g(/-0) 

+ cos ^7r — — /i^ sin (2mq(fi)) sin — /i 

for |/i — 4 1 <C 1. This can be simplified even more to 

F({jl) ~ sin ^7r — \J ^ — cos (2mq(fi)) cos (q(fJ>)) (3-5) 



+ cos (n (tt - V 4 - ) sin (2mg(/i)) ( 2 V 4 ~ M 



where we used the series 

sin Q = ^ V 4 ^ + <?(l 4 - M| 3/2 ) 

with |/i - 4| < 1. Define \i x = 4 - ^ (Z - n) 2 , then 

«0±sin(2m( ? ( W )) QV 4 "^) 
^0 

for Z « n. The eigenvalues of the discretized Laplacian with Dirichlet boundary conditions, defined on the 
real domain [0, 1], are Xi = 4sin 2 (^) with 1 < Z < n — 1. For Z ~ n we have 

A l = 4-^((-nf + 0((l-n) 3 ) 
= M i + 0(a-n) 3 ). 
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So we found that the eigenvalues of —h 2 Lh in the neighborhood of 4 can be approximated by the eigenvalues 
of the Laplacian on the real part of the domain. In the same way we can show that the eigenvalues in 
the neighborhood of \ can be approximated by the eigenvalues of the Laplacian defined on the complex 
contour, by approximating 

«(a0 = \ - + o{\a - 7 Vl 3/2 ). 

Now we are looking for the smallest eigenvalues of —h 2 Lh- Using Taylor series we get 

p( M ) = ^ + 0(M 3 / 2 ) 
<Z(M)=7^ + 0(H 3/2 ) 

so F(p) can be approximated by 

F(p) « sin (n<y/jjL) cos (mjy/Ji) cos 
+ cos (nyfji) sin (777,7 y^) cos 

for p <C 1. We write 7 = 1 + is with < e < 1. This is true for ECS with an angle < 6 1 < J. Then the 
function can be simplified even more to 

F(p) w sin (nyTi + mi^fp) cos(^) (3.6) 

where we used the series 

cos ( 7 ^)=cos(^)+0(M) 

The eigenvalues Xi = ^ ^ of the scaled continuous operator —h 2 L, with Z £ N, are roots of the simplified 
function (3). So the smallest eigenvalues of (3.2) can be approximated by Xi <C 1. □ 

For the Laplace problem (3.2) the spectrum has a typical pitchfork shape. There is a clear complex branch as- 
sociated to eigenvectors located on the complex contour, and a branch closer to the real axis that corresponds 
to eigenvectors located on the real domain. The smallest eigenvalues in the tail of the pitchfork belong to 
the smoothest eigenvectors spread over the entire ECS domain. They lie close to the smallest eigenvalues of 
the continuous ECS operator — L (Figure 5). Indeed, define the complex mesh width h a = R z /(n + m + 1), 
belonging to a straight complex grid connecting and R Z1 and a = h a /h. Then we conjecture for the 
discretized Laplacian —Lh in (3.2) 

a(-L h ) C (6UT) 

where (5 is a strip around the complex line 4p/(ah) 2 (0 < p < po < 1) and % is the interior of the triangle 
P1P2P3 C C, with pi = 4po/(ah) 2 , P2 = 4/h 2 and ps = 4/(jh) 2 . We liberally use the terms, pitchfork 
and tail of the pitchfork to represent the triangular region T and the line segment 6 respectively. For the 
Helmholtz operator with a constant wave number k the pitchfork is shifted in the negative real direction 
over a distance k 2 . 





Remark 2 (Spectrum of smoother ECS transformations). In our analysis we have focused on problems on 
a domain that is complex stretched by the ECS transformation in (2.5). This leads, in ID, to a tridiagonal 
matrix which is constant in the interior and constant in the boundary layer. Both regions are connected by a 
single condition which is the finite difference approximation of the equation at the turning point. For other 
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Figure 5: An illustration of the result in Lemma 3.3. The eigenvalues of the ECS Laplacian discretization 
matrix (•) lie along a pitchfork shape figure, close to the eigenvalues of the same Laplace problem restricted 
to the interior real domain (<) and the complex contour (>) respectively. In the detail view of the area 
around the origin we observe that the smallest eigenvalues practically agree with the smallest eigenvalues of 
the Laplace problem defined on the complex line [0,R Z ] (A). The result is a pitchfork shape: the smallest 
eigenvalues are aligned until they split up into two branches in the point fi\ = 4po/(ah) 2 , with limiting 
points H2 = 4:/h 2 (+) and /x 3 = 4/(jh) 2 (x). (color online) 

complex stretching transformations like quadratic or polynomial scaling (see e.g. Figure 6) the discretization 
matrix is no longer constant and it is much harder to derive theoretical results for the eigenvalues. Numerical 
experiments, however, show a very similar eigenvalue spectrum with a pitchfork. Again, here are some 
numerical eigenvalues, corresponding to smooth modes, that approximated the analytical result, (j7r/R z ) 2 . 
At the pitchfork, the spectrum breaks again into two branches. One branch belongs to eigenmodes that are 
mainly located in the boundary layer and these modes have eigenvalues with a large imaginary part. The 
other branch corresponds to states that are located on the real part of the grid and the eigenvalues will lie 
close to the real axis. 

4 The Model problems 

The actual discrete problems that we solve in the section with numerical experiments, Section 6, are derived 
from three model problems that are representative for break-up problems as they appear in physical systems. 
The discrete formulation is constructed both with the first order Sommerfeld radiation boundary conditions 
as well as with the ECS layers. Therefore, the boundary conditions are not part of the nomenclature. E.g., 
MP1 refers to Model Problem 1 and does not take into account the boundary conditions, which would be 
explicitly mentioned. Likewise MP2 and MP3. Collectively, these model problems are given by the Helmholtz 
equation, 



and are distinguished by the concrete form of the space dependent wave number </>, the right hand side x> 
and the domain size r. 

For a Helmholtz equation on a unit square domain with a wave number cj){x, y) = k 2 , an accuracy condition 
that guards against phase errors polluting the computations [19, 20], requires bounding k 3 h 2 by a small 



{A + <j>(x, y)}u(x, y) = xO, y)\ 



(x,y) e (0,r) 2 



(4.1) 
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Figure 6: A one-dimensional domain with smoother ECS transition in order to preserve the order of dis- 
cretization. The turning point is magnified. 



constant. A similar but less strict constraint that ensures using at least 10 points per wavelength of the 
solution translates to kh < 0.625 [21]. Since we plan to observe solely iterative behavior, we will stick to the 
latter relaxed condition in our experiments. 



4.1 Model Problem 1 (MP1) 

MP1 is the same model problem that forms the basis for the results that appeared in [22]. It is characterized 
by a point source in the center of the domain embodied by a Dirac delta right hand side. In the discrete 
version of the problem, the right hand side is non-zero (= 1) for only one computational node in the scheme, 
which lies at the center of the domain. The wave number is constant in MP1, and the domain is a unit 
square. 

Specification of MP1: (4.2) 
4>(x,y) = k 2 , k G M (constant wave number) 

x(x, y) = S(x, y) 5(x, y) (Dirac's delta function) 

r = 1. 

A plot of MP1 with both type of boundary treatments is given in Figure 7. 



4.2 Model Problems 2 and 3, (MP2, MP3) 

MP2 and MP3 are Helmholtz model problems with strongly varying wave numbers, and therefore pose a 
tough benchmark for an iterative approach. These problems originate from large scale Helmholtz problems 
that appear in the simulation of Schrodinger's equation for single and multiple ionization of atoms and 
molecules [4]. The dynamics of two positively charged electrons i*i and r 2 in the field of the negatively 
charged nuclei need to be modeled and solved. Usually the electrons move, due to their mass difference, 
much faster than the nuclei and the latter are usually taken fixed in space. This leads to a Schrodinger 
equation with a six-dimensional wave function ^(ri,^), which is often expanded in spherical coordinates 
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(a) MP1 with ECS layers (<9 7 = f ) 



(b) MP1 with Sommerfeld BC 



Figure 7: The real part of the solution of MP1 with k = 160, discretized on a grid having 256 interior points 
per dimension. This suffices to meet the accuracy condition, (color online) 

about the center of the molecule with the z axis along the axis of the molecule. The expansion is 



where (pi, Oi) and (p2 : 2 ) are the spherical coordinates of the first and the second electron and ^ denotes 
the two angles in the spherical coordinates. This expansion leads to a very large number of coupled 2D 
problems, where ^i 1 m lj i 2 m 2 (pi, P2) is the solution of the 2D problem for particular integer values of Zi, mi, 
I2 and 777,2. It describes the wave as a function of the distances p\ and p2 of both electrons to the center of 
the molecule. The coupled equation has a block structure and the differential operators only appear in the 
diagonal blocks since Yi m (Q) are eigenfunctions of the angular part of the Laplacian operator in spherical 
coordinates. 

In the work [4], the resulting linear systems for a molecule are solved iteratively. The problem was pre- 
conditioned by inverting the diagonal blocks with the direct sparse solver SuperLU, which is based on the 
left-looking supernodal method [23] and was employed on a massively parallel computer. The current work 
aims to replace this direct solver with an iterative alternate based on preconditioning. The main motivation 
for studying the numerical properties of the solver for the model problems is that the approach in [4] can not 
be used for solving systems with three or more particles. The diagonal blocks that need to be inverted in 
this extended case are at least 3D problems, and the resulting storage and computational complexity grows 
out of reach for the current computational infrastructures. 




("2), 



Specification of MP2: 



(4.3) 




< k < 5, 



< v < 10 



x(x,y) = 



1 



e x 2 +y 2 



r = 50. 



Specification of MP3: 



(4.4) 



14 



<j>(x, 2/) = - + - + k 2 , <k <5 

x y 

50 <r < 200 

These model problems are representative of the 2D problem that appear when l\ = and I2 = 0. The 
coordinates x and y should be interpreted as radial variables p\ and p<i. 

Homogeneous Dirichlet boundary conditions stay fixed at the south and the west edges of the domain where 
pi and p2 are zero. On the east and the north edges where p\ or p2 are large, absorbing boundary conditions 
have to be used. We therefore toggle between the first order Sommerfeld BC and the ECS layers on these 
two edges, and provide numerical results with both. For quality calculations that reproduce the physical 
experiments, higher order absorbing boundary conditions need to be used, but for the purposes of the 
current paper these low order boundary conditions are sufficient. The boundary conditions that MP2 and 
MP3 employ are given by Equation (4.5). 

u (Q->y) — u(x,0) = Homogeneous Dirichlet BC, south/west edges (4.5) 
— = — iku Sommerfeld BC, east/north edges 

(4.6) 

= ECS layers, east/north edges 

For v = 7 and k — 2, a plot of the solution of MP2 appears in Figure 8. In this particular case, the minimum 
grid size required for an acceptable resolution of the solution is 341 2 , closest to which the most convenient 
practical grid size is 384 2 from a multigrid perspective. Later in Section 6.2, we describe how we evaluate 
the minimum grid size for MP2 and MP3. The solution has evanescent waves for values of v > 2.73, which 
are damped exponentially on these edges. These evanescent waves correspond to single ionization break-up 
reactions [3]. 

In MP3, <j)(x, y) has a singularity at the origin. Unlike MP2, the solution of MP3 always has evanescent 
waves at the south/ west edges, regardless of the choice of the parameters r and k. For r = 90, and k = 2, 
the minimum interior grid size required is around 1024 2 and the solution is depicted in Figure 9. 




Figure 8: The real part of the solution of MP2, with v = 7 and k = 2. For values of v > 2.73, evanescent 
waves form near Dirichlet edges, (color online) 

This completes the description of the model problems, which we solve iteratively in Section 6.2. 
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(a) MP3 with ECS layers (0 7 = f ) 



(b) MP3 with Sommerfeld BC 



Figure 9: The real part of the solution of MP3, with r = 90 and k = 2. Evanescent waves are visible near 
Dirichlet boundaries. The wavelength is quite small and implies a huge grid size to comply with accuracy 
requirements, (color online) 

5 Multigrid and preconditioning 

In this section we discuss the effect of the spectral properties of the indefinite Helmholtz operator on the 
convergence of iterative processes such as multigrid. The spectrum Xn h G cr(Hh) of the discrete Helmholtz 
operator and that of the discrete Laplacian on an ECS grid, i.e., \l h £ &(Lh)i varies only up to a real 
constant — /c 2 , which defines the distance by which an eigenvalue shifts westwards to render A# h . We 
saw in Section 3 that Xn h is constrained to an area consisting of a straight line segment, and a region bounded 
by a triangle, the so-called pitchfork. We observe that for small values of hk, the smallest eigenvalues lie on 
the tail of the pitchfork. The eigenmodes corresponding to these small eigenvalues are the standing waves 
that cover both the real and the complex part of the domain. It is important to note that none of the small 
eigenvalues can be exactly equal to zero because we have shown that each of the corresponding modes must 
possess at least a non-zero imaginary part. Depending on the magnitude of the real wave number, this can 
potentially lead to a very large (but bounded) condition number and thus confirms that the discrete problem 
is ill-conditioned. 

5.1 Multigrid 

In this section, we assume familiarity with basic geometric multigrid. See [24, 25, 26] for a quick access. 
Here, we briefly skim through the multigrid difficulties in solving an indefinite Helmholtz problem. 

The first aspect of multigrid that requires attention in the context of indefinite linear systems, is the absence 
of a pointwise smoothing procedure. For a given discrete operator Mh, a strict condition on the so-called h- 
ellipticity measure Eh(Mh) [27, 26], viz, Eh(Mh) > formally implies the existence of a pointwise smoothing 
process. Circumventing the details, it suffices to mention here that the /i-ellipticity of the discretized indefinite 
Helmholtz operator is very close to zero for interesting values of the wave number, and therefore, common 
stationary methods do not amply relax the error to be representable on the coarse grid. 

The other troublesome multigrid aspect that merits attention is coarse grid correction. To see this, imagine 
that we increase the mesh width while keeping k constant. The rightmost eigenvalues of the discrete operator 
on the finest level lies near A/h 2 — k 2 and cannot come near zero without strictly violating the posted accuracy 
condition hk < 0.625. However, this is more likely to happen within a multigrid cycle where the same operator 
is re-discretized on the coarser grids, each with a larger mesh width. This effect can lead to resonant behavior 
on a level where A/h 2 — k 2 « 0. This leads to a severe degradation of multigrid performance. This issue is 
also well-known and discussed in papers, such as [28, 21]. 
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An important point to observe in this context, is that smoothing is much less of a trouble than coarse grid 
correction in the case of a multigrid solution of an indefinite problem. For moderate wave numbers in MP1 
such as k = 40, we see that multigrid still converges after a careful choice of components. For example, 
we conducted a test on this problem with k = 40 and interior grid size TV = 64 in both dimensions; with 
ECS layers on all edges. The components were ILU(O) with u = 0.3, F(l,2)-cycle, FW-restriction, bilinear 
prolongation and Galerkin formulation of the coarse grid operator. For this problem multigrid converged 
in 21 iterations. However, we also saw that for twice the grid size and twice the wave number, the same 
algorithm failed to converge with any combination of components. When we introduced a slight damping in 
the wave number, k = (80 — 0.05z), multigrid converged again in 35 F(l,l)-cycles, and thus confirms that it 
might be used for approximate preconditioner solves in a Krylov setup. Note that these choices of k and N 
are not ideal from an accuracy perspective, however, the above example is useful to understand the multigrid 
performance. 

5.2 A short overview of the complex shifted Laplacian (CSL) preconditioner 

The idea of preconditioning the Helmholtz problem with its (slightly) damped version as published by 
Erlangga et al. in [13, 21, 22] is founded on avoiding the diverging behavior of multigrid for the original 
indefinite problem. The damping is brought about by a complex shift of the Laplacian, which evidently shifts 
its spectrum away from the origin. As a result, the discrete problems formulated with the shifted Laplacian 
can be tackled by multigrid. This preconditioning is perfectly applicable in the ECS context as well. The 
continuous version of the preconditioner for the Helmholtz problem on the ECS domain [0, 1] U C C 

is given by 

M CSL = -(£ + p 2 k 2 ) in [0, 1] U [1, R z ] C C, (5.1) 

with a complex shift (3 2 = e\ + is<2 G C. It is also important to know that in the comparison with [21], the 
complex shift j3\ — ifa in [21] is equivalent to e\ + is^ in this paper. 

Since M CSL = — (L + f3 2 k 2 ) and H — —(L + k 2 ) share the same eigenvectors it is easy to see that the 
eigenvalues of the continuous preconditioned system (M csx ) -1 i7 lie on a circle. Indeed, the spectrum is 
given by the linear fractional transformation LF(fi) = J^2 k 2 that maps the complex line of eigenvalues Xl 
of the Laplacian to the circle through 1 ^2^2 fc 2 and 1. 

The linear fractional transformation maps the negative imaginary half plane to the interior of the complex 
circle through the complex values 0, 1 and 1//3 2 . So the spectrum of the discretized operator is mapped 
to this region. More specifically, each line of the triangle that bounds the eigenvalues is mapped onto a 
segment of a circle that lies inside this region. As a result, the eigenvalues of the preconditioned discrete 
system (M^ SL )~ 1 Hh lie away from the origin, inside a banana shaped figure that is the image of the different 
branch lines in the pitchfork and the tail from Figure 5 as illustrated in Figure 10. The preconditioned system 
can be solved much more conveniently with Krylov subspace solvers by employing very few multigrid cycles 
for approximate preconditioner inversion. 

5.3 The complex stretched grid (CSG) preconditioner 

We readily see that the indefinite spectrum may be shifted favorably by an alternate strategy, which is 
more focused towards the ECS formulation. Instead of scaling the wave number, we keep it unchanged and 
scale the discretization grid. To see this, imagine discretizing the one-dimensional complex shifted Laplacian 

on a real interval with constant mesh width h G M, and note that 
(f3 2 k 2 )h 2 = k 2 (f3 2 h 2 ). The left term in this last equality appears in the discretization of M CSL , while the 
right term can be interpreted as a quantity that appears in the discretization of the Helmholtz problem 
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Figure 10: The pitchfork of the original Helmholtz operator in Figure 5 is mapped to a banana (dotted lines), 
with the eigenvalues (•) inside. The lines corresponding to the eigenvalue problem on the real domain and 
the complex contour are mapped to a circular arc (<) through and a smaller circular arc (>) respectively. 
They enclose the preconditioned spectrum and the image of the intermediate complex line associated to the 
eigenvalue problem on the complex line [0, R z ] (A), (color online) 



Mcsg — ~j^2 — k 2 defined on a straight complex line with constant mesh width f3h G C. Indeed 

Mf SL u h = ~(^L h + (3 2 k 2 )u h = b h & M% SG u h = -(p^L h + k 2 )u h = ^b h , (5.2) 

so the system M GSL Uh = bh yields the same solution as M GSG Uh — bh/fi 2 . In this way we found the 
equivalent complex stretched grid (CSG) preconditioner M CSG . In this context we will denote the angle of 
P in the complex plane as Op. 

The same argument still holds on ECS domains where the contour mesh widths are already complex and 
for inhomogeneous wave numbers as present in MP2 and MP3. This approach offers extra possibilities to 
explore. Instead of scaling, for example, the entire spectrum away from the origin, only the problematic 
branch of eigenvalues close to the real axis can be scaled deeper into the complex plane. In other words, 
only the interior part of the grid is scaled with /?, while the complex contour stays the same. In general, we 
can build a preconditioner by defining the original Helmholtz equation on a convenient domain, given by a 
coordinate transformation 

J x + ifa(x), x G 
z[x) = < 

I x + ifr(x), x G T, 

with fn, fr G C 2 increasing (e.g. linear, quadratic, . . . ) and lim fn(x) = lim fr(x) ^ 0. 
Applied to (2.5) we get 

M CSG u(z) = -^u(z) in [0, z(l)} U [z(l), z(R)] C C, (5.3) 

with homogeneous Dirichlet boundary conditions in z(0) = and z(R) = R z . Independent of the complex 
contour T z = we now have the interior region Q z = [0,2(1)] complex as well. 

In Figure 11 the interior region is chosen along the line connecting the two boundaries and R z of the 
original ECS domain. As a result the interior mesh width is scaled from h to /3ft,; the complex contour has no 
extra scaling, the mesh width 7ft is preserved. Figure 13 displays the resulting preconditioned system with 
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CSG (right) versus the CSL (left) with the same scaling in the wave number, from k to (3k. CSG leads to a 
similar C-shaped spectrum, away from the origin, favorable for Krylov methods. However, experiments show 
that the spectrum of the preconditioner (see Figure 12) is still bad for our current multigrid configuration. 
The numerical experiments with the CSG preconditioning method in the next section all use grids that are 
equally scaled over the entire domain, i.e. the interior part and complex contour. The CSG preconditioner 
with scaling factor (3 is then related to the equivalent CSL preconditioner with scaled wave number /3k, this 
means a complex shift (3 2 , as in equation (5.2). 
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Figure 11: The CSG domain. The domain for M CSG (line) is complex stretched in the region of interest 
[0, 1] as well, such that it aligns with the line connecting the two boundaries and R z of the original ECS 
domain (dashed line). 




Figure 12: The eigenvalues of the discrete preconditioner M GSG on the domain in Figure 11. The eigenvalues 
lie in a narrow pitchfork (•) that lies in the bottom half of the pitchfork of the original matrix — (dotted 
lines) from Figure 5. By choosing the domain as in Figure 11, the top branch of the new pitchfork lies along 
the middle line associated to the eigenvalue problem on the complex line [0, R z ], i.e. the line of eigenvalues 
of the continuous problem. 



6 Numerical experiments 

In this section we supplement the theoretical development with numerical experiments carried out on the 
model problems with both kinds of boundary treatment as discussed earlier. We provide a brief comparison 
of different multigrid components, which allowed us to choose the best set for the approximate multigrid 
inversion of the preconditioning operator during a Krylov solve. The details of the solver and results from the 
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(a) Complex shifted Laplacian 



(b) Complex stretched grid 



Figure 13: Left: The banana shaped spectrum of the CSL-preconditioned Helmholtz system (M^ 75 ' L ) _1 i^/ l 
in Figure 10 grows slightly towards the origin for increasing k. The same effect was observed for Som- 
merfeld radiation conditions by Erlangga et al. in [22]. Right: The eigenvalues of the equivalent discrete 
preconditioned Helmholtz system (M^ SG )~ 1 Hh with the complex stretched grid preconditioner defined on 
the domain in Figure 11, lie on a C-shaped figure, away from the origin, with a similar dependence on k. 

numerical experiments are displayed in tables and figures for easy access. For the numerical experiments we 
use the Shortley-Weller finite difference discretization (3) applied to a cell-centered mesh topology. The cell- 
centered mesh topology is chosen because multigrid is slightly more convenient with this choice for general 
Robin-type first order derivative boundary conditions (henceforth BC), of which the Sommerfeld radiation 
BC are a special case. Moreover, the stencil is the same as for the vertex-centered case, so the results are 
easily carried over. 

6.1 Choice of Multigrid Components 

To select a set of multigrid components from the different available choices, we take MP1 with k = 80, 
and use a grid size of 128 2 obeying the accuracy constraints. ECS layers are used, and defined by scaling 
the mesh widths in these layers by e %e ^ with the angle # 7 = 7r/6. In this case, most of the spectrum lies 
in the 4 th quadrant of the complex plane, i.e., except the eigenvalues responsible for making the linear 
system indefinite. Through a negative imaginary shift of the Helmholtz operater equal to —0.2, we push the 
spectrum adequately towards the 4 th quadrant, thus transforming the linear system so that it is now nearly 
negative semi- definite. It is imperative to comprehend that with Sommerfeld BC, the major part of the 
spectrum of the original problem is in the 1 st quadrant, and therefore with Sommerfeld BC, a preconditioner 
formed through a positive imaginary shift will make sense. 

From Figure 14, we immediately recognize the set of multigrid components that work best in the present 
context. ILU(O) smoothing, Full Weighted (FW) averaging as restriction, bilinear interpolation as prolon- 
gation, and the Galerkin formulation for the coarse grid operator. These components are harnessed within 
a V(0,l)-cycle and yield an algorithm that does an excellent job on the preconditioner. Therefore, these are 
the components that we will invariably use in Section 6.2. For the sake of completeness, we also checked 
out multigrid performance with various combinations of cj-Jacobi, with an under-relaxation of 0.5, employed 
in F-cycles. The comparison was very thorough, and included tests with the Four-Point averaging [26] as 
restriction, as well as with direct discretization on the coarse grids. The results are presented in Figure 14, 
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' co-Jacobi/0.5/F(1,2)/FW/GCG 
u w-Jacobi / 0.5 / F(1 ,2) / FP / DCG 
— ILU(O) / 1 .0 / V(0,1 ) / FW / GCG 
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Number of multigrid cycles 



(a) Defect reduction / mg cycles 




CPU time (seconds) 
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Figure 14: Multigrid performance on the preconditioner formed through the CSL technique (f3 = (1, —0.2)) 
for MP1, with k = 80, 128 interior points plus 2 x 32 points in the ECS layers along each dimension (on 
either sides). In all the cases, prolongation was done through bilinear interpolation. The legend incorporates 
the smoother, the relaxation parameter, the restriction method, the cycle type, and the coarse grid operator. 
DCG stands for Direct Coarse Grid operator, while GCG stands for the Galerkin Coarse Grid operator. 



k 

V 


1 


2 


3 


4 


5 





0.625 


0.312 


0.208 


0.156 


0.125 


1 


0.360 


0.255 


0.188 


0.147 


0.120 


2 


0.279 


0.221 


0.173 


0.140 


0.116 


3 


0.236 


0.197 


0.161 


0.133 


0.112 


4 


0.208 


0.180 


0.151 


0.127 


0.109 


5 


0.188 


0.167 


0.143 


0.122 


0.105 


6 


0.173 


0.156 


0.136 


0.118 


0.102 


7 


0.161 


0.147 


0.130 


0.114 


0.100 


8 


0.151 


0.139 


0.125 


0.110 


0.097 


9 


0.143 


0.133 


0.120 


0.107 


0.095 


10 


0.136 


0.127 


0.116 


0.104 


0.093 



Table 1: Maximum mesh width limits for MP2 



are self explanatory, and indicate the most viable set of components very clearly. 



6.2 Numerical Experiments 

In this section before putting up the results of the numerical experiments, we first sort out the minimum 
grid size requirement for MP2 and MP3, which are Helmholtz problems with strongly varying wave numbers. 
As the wave numbers are infinite at the origin, we take into consideration the highest discrete wave number 
that the discretized problem can attain under the cell-centered mesh topology. For reasons of brevity, the 
complete analysis is not shown here. We however, brief the steps involved in the analysis. First, we transform 
the problem into a new coordinate system; so that x = x/r, and y = y/r. This new system is in the unit 
square domain and is therefore dimensionless. The next step is to evaluate the supremum of the transformed 
wave number, which occurs at the origin for the model problems. For MP2, this is straightforward. We use 
this value in the accuracy condition and are led to the maximum mesh widths that must be obeyed. For 
MP2 this is given in Table 1. 

Here, we will solve two test-cases of MP2. One is characterized by v = 7 and k = 2, while the other is 
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k 

r 


1 


2 


3 


4 


5 


V 50 < r < 200 


0.095 


0.089 


0.082 


0.075 


0.068 



Table 2: Upper limit on the mesh width for MP3 with respect to k 









MG perform., 


MGP 


MGP 


MGP 




BC 


Preconditioning operator 


on precond. 


Bi-CGSTAB 


IDR(4) 


IDR(8) 




Conv., # iter 


matvec 


matvec 


matvec 








cputime 


cputime 


cputime 


cputime 




O 
PQ 


CSL, with shift = -1 + 0.2i 


0.31, 12 
0.52 sec. 


79 
4.75 sec. 


72 
4 sec. 


76 
5.35 sec. 


o 


Lmerfel 


7T 

CSG, with angle = — — 
28 


0.24, 10 
0.47 sec. 


74 
4 sec. 


70 
3.98 sec. 


72 
4.40 sec. 


to 

1—1 

II 


a 

o 

CO 


CSL, with shift = -e~ 2l7T/28 = -0.97 + 0.22z 


0.26, 11 
0.48 sec. 


74 
4.78 sec. 


69 
4.26 sec. 


72 
5.20 sec. 


MPl: 


ECS Layers 


CSL, with shift = -1-0.22 


0.29, 12 
1.13 sec. 


59 
9.15 sec. 


62 
9.11 sec. 


60 
9.61 sec. 




7T 

CSG, with angle — — 
28 


0.24, 10 
1.12 sec. 


58 
9 sec. 


62 
9.23 sec. 


56 
9.62 sec. 




CSL, with shift = - e 2tn/28 = -0.97 - 0.22z 


0.24, 10 
1.10 sec. 


58 
9.04 sec. 


59 
9.06 sec. 


54 
9.30 sec. 



Table 3: Experimental results - Iterative solution of MPl. One V(0, 1) multigrid cycle is used for approximate 
inversion of the preconditioners. This preconditioning is used with Bi-CGSTAB, IDR(4) and IDR(8). The 
interior domain is a unit square, discretized with 256 points along both dimensions. With ECS layers, there 
are 2 x 64 additional points per dimension, belonging to the layers. With ECS formulation each edge of the 
domain is endowed with an absorbing layer. 



characterized by v = 1 and k = 4. For both of these test-cases, the minimum grid size (from a multigrid 
perspective) is 384 2 . 

For MP3, we immediately see that the origin cannot be substituted into the dimensionless wave number (due 
to the singularity). We observe that (h x /2,h y /2) gives the largest wave number that the discrete problem 
can attain. Using this in the accuracy condition results in a quadratic equation, which we solve to get the 
mesh sizes. These are scaled back to (0,r) 2 and are given in Table 2. 

We form two test-cases with MP3 as well. The first one is characterized by r = 90 and k = 2 and requires a 
minimum interior grid size of 1024 2 . The second test case is much more severe. It is formed by r = 150 and 
k = 4, and requires a minimum interior grid size of 2048 2 for an acceptable resolution. 

Remark 3 (Reading and comparing the experimental results). The experimental results are summarized 
in Tables 6.2, 6.2, and 6.2. Table 6.2 accounts for one, and Tables 6.2 and 6.2 account for two problems 
(each) derived from using different values of the parameters in the model problems. This is mentioned 
vertically in the first column of each table. Each of these problems give two different discrete versions, when 
formulated once with the Sommerfeld approximation (on the boundaries), and second, when formulated 
with ECS layers. This is also vertically marked in the tables. Next, each discrete formulation is solved with 
Krylov methods, using three preconditioners; (1) the CSL preconditioner with the shift k 2 — » f3 2 k 2 where 
ft 2 = E\ + iS2 = 1 + £2 (with 62 having the smallest absolute value for which the chosen multigrid method 
converges), (2) the CSG preconditioner having the scaling h — » j3h where f3 = e z6 P (with 0p the smallest 
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angle for which the chosen multigrid method converges), and finally, (3) the CSL preconditioner again, with 
a complex shift equal to j3 2 = e 2z6>/3 (with Op as used with the CSG preconditioner). Each row in the tables 
accounts for one of the above mentioned preconditioners. Each row (after specifying the preconditioner) lists 
the following information: 

1. Multigrid performance on the preconditioner taken as a standalone problem. Com;., is the average 
convergence factor per cycle, and #iter is the number of cycles consumed to reduce the relative 
residual by 7 orders of magnitude. 

2. The number of multigrid preconditioned (MGP) matvecs employed by Bi-CGSTAB. Note that each 
iteration of Bi-CGSTAB consists of two such matvecs. We chose matvecs over conventional iterations 
to have a better idea of the total number of multigrid cycles used, as well as to have a fair comparison 
with IDR(s). 

3. The number of multigrid preconditioned matvecs employed by IDR(4). 

4. The number of multigrid preconditioned matvecs employed by IDR(8). 

Within each row, the three Krylov methods can be compared against one another with the same precondi- 
tioner. Within the stacks containing three rows each, the relative performance of the three preconditioners 
on an identical discrete problem can be checked. 

In this paper we use one multigrid V(0, l)-cycle for each preconditioning step that is involved in the solver. 
Convergence of the Krylov solver is determined by a check on the relative residual going below a tolerance 
value of 10 -6 , i.e., the algorithm stops after the m th iteration if: 

w < 10 (fL1) 

|| 2 is the defect (measured in the discrete I/2-norm) after the i th solver iteration. 

It is a well-known fact that due to indefiniteness multigrid cannot work directly with the Helmholtz problem 
as a solver. The CSL and the CSG preconditioners are attempts to maneuver the spectrum slightly so that its 
indefiniteness can be reduced, and that the eigenvalues may be slightly shifted so as to bring them closer to 
the positive or to the negative definite regimes in the complex plane. Roughly speaking, the CSL translates 
the spectrum, while the CSG rotates it to accomplish the objective. A preconditioner, therefore, can also be 
build as a hybrid between the CSL and the CSG approaches, i.e., with a general CSL shift e\ + 1S2 combined 
with a general CSG rotation angle Op. However, during experimentation, the latter preconditioner did not 
prove any better than either of these two approaches used in isolation, and is hence not presented. 



Remark 4 (On choosing the CSL and the CSG parameters). An automation can be set up which starts 
from a given imaginary shift for the CSL preconditioner, or a given angle for the CSG precondtioner, and 
monitors the residual norm obtained after successive multigrid cycles. In such an automatic routine, absolute 
values of the shift size or the angles may be reduced to a benchmark for which the given multigrid method 
just converges (say in 10-20 cycles). Parametrized with this shift (or angle), the preconditioner may be used 
in the Krylov method with approximate solves for preconditioning. Note that such a selection routine may 
only run once, and decide upon the shifts to be employed for all later Krylov iterations. In this paper, 
however, we just resort to doing the above described process manually. 

We first observe from Table 6.2 that the performance of the CSL preconditioner with both shifting strategies 
is similar to the performance of the CSG preconditioner. The multigrid inversion of the preconditioners as 
well as the overall numerical solution method are very good for MP1. 
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Table 4: Experimental results - Iterative solution of two different discrete problems obtained from MP2 
(by using different values of v and k). One V(0, 1) multigrid cycle is used for approximate inversion of the 
preconditioners. This preconditioning is used with Bi-CGSTAB, IDR(4) and IDR(8). For both problems, the 
interior domain is a square of 50 units, discretized with 384 points along both dimensions. With ECS layers, 
there are 128 additional points per dimension, belonging to the layers. Contrary to MP1 the formulation 
with ECS Layers only has these layers on two edges of the domain, the north and the east. 

In Table 6.2, we have solved two test-cases of MP2 with different characterizing parameters. The values of 
v and k distinguish the test-cases. From Table 1, we read that the mesh width requirement for both these 
test-cases is the same (0.147), and therefore they can be solved on an identical grid of interior size 384 2 . We 
observe that although the supremum of the wave number in the domain is identical for both the test-cases, 
the second one takes twice the time needed to compute the solution of the first one, for MP2. 

Remark 5 (Smooth ECS transition). Plausibly, the sharp rotation of the linear ECS contour may not 
be very desirable for the discretization of some applications. For this reason, we checked out the iterative 
performance of the preconditioners and the solvers, for MP2, formulated with an ECS contour that rotates 
gradually in 256 very small equally sized angles. The one-dimensional analog of the domain is shown in 
Figure 6. The performance of multigrid (on such preconditioners) as well as the Krylov methods turned out 
very similar to that listed in the tables. 
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Table 5: Experimental results - Iterative solution of two different discrete problems obtained from MP3 (by 
using different domain sizes and values of k). One V(0, 1) multigrid cycle is used for approximate inversion 
of the preconditioners. This preconditioning is used with Bi-CGSTAB, IDR(4) and IDR(8). The interior 
domain is a square of r units. For r = 90, k = 2 (upper six rows), the problem is discretized with 1024 
interior points along both dimensions. In these problem, when ECS layers are used, there are 256 additional 
points per dimension, belonging to the layers. For r = 150, k = 4, the problem is discretized with 2048 
interior points per dimension. In the ECS formulation there are 512 extra points (per dimension) in these 
layers. For both problems, ECS layers are only required, and used, at the north and the east edges of the 
domain. 

The results from experiments on the harder model problem, i.e., MP3 whose spectrum is more indefinite 
compared to the other model problems here, are laid out in Table 6.2. We clearly see an advantage in the 
number of matvecs on problems with the ECS layers against the Sommerfeld BC. However, each matvec of 
an ECS problem takes longer due to the additional grid points in the layer, hence this advantage is not seen 
in the CPU time. Figure 15, details the convergence history of the first test-case of MP3. This is important 
in order to check out any possible stagnation trend in the convergence. We find however, that in all the cases 
Bi-CGSTAB and IDR(4) show a well-matched convergence behavior. IDR(4), however, seems to be slightly 
faster in comparison with Bi-CGSTAB as it avoids the slight initial stagnation evident with Bi-CGSTAB in 
the figure. With a comparison between two Krylov methods, it is often also interesting to observe both the 
CPU time as well as the solver matvecs due to possible differences in the subspace minimization strategies 
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(b) IDR(4) on MP3 with Sommerfeld BC 
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(c) Bi-CGSTAB on MP3 with ECS layers 
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(d) IDR(4) on MP3 with ECS layers 



Figure 15: Bi-CGSTAB and IDR(4) defect reduction history with the proposed CSG preconditioner as well 
as the CSL preconditioner for MP3 with 2048 2 points in the interior of the domain. ECS layers (where used) 
have added 512 points per dimension in the exterior layers. Expected convergence traits for Bi-CGSTAB 
and IDR(4) are clearly visible, (color online) 



in different methods. This reveals that although for many test-runs IDR(8) reports lesser matvecs than 
IDR(4), it really does not provide any concrete enhancement. Evidently, the reduction in matvecs is easily 
offset by the increase in CPU time (plus an associated increase in storage which is not shown). 

We also checked out the solver performance for a larger MP3 problem with size 4096 2 (not reported in the 
tables). The CPU time already runs over 4 hours for a problem of this size, although the performance of 
the iterative method is of the same quality as the other problems. This is due to the huge complexity of the 
test. 



Remark 6 (Testing the situation with GMRES). We ran multigrid preconditioned GMRES on MP2, with 
v = 1,/c = 4 (ECS formulation). Multigrid was used to approximately invert the CSG preconditioning 
operator exactly as specified in the second last row of Table 6.2, during each GMRES iteration. GMRES 
reported 103 iterations to reduce the relative residual by seven orders of magnitude, and reported the 
consumed CPU time as 175 seconds. The same level of accuracy can be reached with Bi-CGSTAB or IDR(s) 
(as depicted in Table 6.2 in roughly one fifth of the time required for GMRES. GMRES therefore qualifies 
well for analytic purposes, or where the problem size is small. 



Remark 7 (Platform specification and storage scheme). These experiments were performed serially on an 
Intel Xeon (8-core) with 32 gigabytes of RAM. We used Matlab v7 as the testing platform. Some processes 
such as matvec computations, and backslash inversion of some matrices in the multigrid heirarchy were 
performed in parallel using all the 8 processors. We implemented some linear algebra routines in mex-C 
to speed up the computations, and used the Compressed Row Sparse Format for storage of the matrix 
data. Matlab uses Compressed Column Sparse Format for such needs. Both the CSL, as well as the CSG 
preconditioners were stored as sparse matrices in either of these formats. 
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7 Conclusions and outlook 



In this paper, we have studied the Helmholtz problems that arise in mathematical models for single and 
double ionization of atomic and molecular systems. The problems typically have regions in space where the 
wave number can be large and absorbing boundary conditions are often implemented with complex stretched 
grids. 

We developed and analyzed the iterative properties of the exterior complex scaled (ECS) absorbing boundary 
layers for the indefinite Helmholtz equation. We have analyzed the spectral properties of the discrete problem 
formulated with ECS layers and found bounds around the spectrum of the Helmholtz operator for constant 
wave numbers. These bounds were derived for a finite difference Shortley-Weller discretization and linear 
exterior scaling. Although the theoretical estimates are limited to this model, numerical tests suggest that 
they are valid for quite general cases. 

An alternative preconditioner to the complex shifted Laplacian (CSL) is introduced where instead of shift- 
ing the wave number, the grids are given a complex scaling. We call this complex stretched grid (CSG) 
preconditioning. We introduced two new benchmark problems that are derived from break-up problems in 
quantum mechanics and have strongly varying spatially dependent wave numbers. They provide a tough 
benchmark for future development of iterative Helmholtz solvers. 

The CSG and the CSL preconditioners are related and perform similarly for most problems. The precondi- 
tioner inversion is performed approximately by a geometric multigrid method, based on ILU(O) relaxation, 
V(0,l)-cycles, FW restriction, bilinear prolongation, and the Galerkin coarse grid operator. Different numer- 
ical experiments with the CSG preconditioner on the model problems, show that our multigrid method is 
more stable for problems with constant wave numbers. With spatially dependent wave numbers we see that 
our multigrid method requires a higher damping of the operator to render the indefiniteness manageable. 
However, with greater damping, the spectra of the preconditioner and the operator get farther apart and 
this prolongs Krylov convergence. Although not listed, we tried different alternatives for multigrid precon- 
ditioning, such as F-cycles and more than one V-cycle, however, the best results from a CPU time aspect 
are shown. 

In future, we intend to explore two different avenues to improve the solver for the model problems introduced 
here. The main problem is the reduction of CPU time for problems formulated with ECS layers. This can be 
brought about by first reducing the complexity of the problem by applying adaptive mesh refinement (AMR) 
techniques for discretization, say, near the evanescent layers of the solution. And secondly, by identifying 
a way of shifting only the most problematic part of the spectrum, since every cluster of eigenvalues that is 
needlessly shifted can significantly decrease the performance of the preconditioning. This may results in a 
better solver than the current state-of-the-art, and is thus our future goal. 
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